function householder_qr()
% HOUSEHOLDER_QR Householder QR分解演示
% 
% 演示基于Householder反射的QR分解算法

fprintf('=== Householder QR分解 ===\n\n');

%% 1. Householder反射基础
fprintf('1. Householder反射基础\n');
fprintf('Householder矩阵 H = I - 2*v*v^T，其中 ||v|| = 1\n');
fprintf('H是对称、正交、对合的矩阵\n\n');

% 演示Householder反射的性质
v = [1; 2; 3];
v = v / norm(v);  % 归一化
H = eye(3) - 2 * (v * v');

fprintf('单位向量 v = [%.4f; %.4f; %.4f]\n', v);
fprintf('Householder矩阵 H =\n');
disp(H);

% 验证性质
fprintf('H的性质验证:\n');
fprintf('对称性 ||H - H^T||: %e\n', norm(H - H', 'fro'));
fprintf('正交性 ||H*H^T - I||: %e\n', norm(H * H' - eye(3), 'fro'));
fprintf('对合性 ||H^2 - I||: %e\n', norm(H^2 - eye(3), 'fro'));

%% 2. 向量的Householder反射
fprintf('\n2. 向量的Householder反射\n');
fprintf('将向量x反射到与第一个坐标轴平行\n\n');

x = [3; 4; 5];
fprintf('原始向量 x = [%g; %g; %g]\n', x);
fprintf('向量长度: %.6f\n', norm(x));

% 构造Householder向量
[v_house, beta] = householder_vector(x);
fprintf('Householder向量 v = [%.6f; %.6f; %.6f]\n', v_house);
fprintf('参数 β = %.6f\n', beta);

% 应用反射
H_x = eye(3) - beta * (v_house * v_house');
y = H_x * x;
fprintf('反射后向量 y = [%.6f; %.6f; %.6f]\n', y);
fprintf('理论结果: [±%.6f; 0; 0]\n', norm(x));

%% 3. Householder QR分解算法
fprintf('\n3. Householder QR分解算法\n');
fprintf('逐列消除下三角部分，得到上三角矩阵R\n\n');

% 测试矩阵
A = [1, 2, 3;
     4, 5, 6;
     7, 8, 9;
     10, 11, 12];

fprintf('原始矩阵 A =\n');
disp(A);

% Householder QR分解
[Q, R] = householder_qr_decomp(A);

fprintf('Q矩阵 =\n');
disp(Q);
fprintf('R矩阵 =\n');
disp(R);

% 验证分解
QR_product = Q * R;
fprintf('验证 Q*R - A 的误差: %e\n', norm(QR_product - A, 'fro'));

% 验证Q的正交性
QtQ = Q' * Q;
fprintf('Q的正交性 ||Q^T*Q - I||: %e\n', norm(QtQ - eye(size(QtQ)), 'fro'));

%% 4. 紧凑存储形式
fprintf('\n4. 紧凑存储形式\n');
fprintf('将Householder向量存储在R的下三角部分\n\n');

[R_compact, householder_data] = householder_qr_compact(A);

fprintf('紧凑存储的R矩阵 (包含Householder向量):\n');
disp(R_compact);
fprintf('β参数: [%.6f, %.6f, %.6f]\n', householder_data.beta);

% 从紧凑形式重构Q
Q_reconstructed = reconstruct_Q_from_compact(R_compact, householder_data);
fprintf('重构的Q矩阵与原Q的差异: %e\n', norm(Q_reconstructed - Q, 'fro'));

%% 5. 与Gram-Schmidt的比较
fprintf('\n5. 与Gram-Schmidt的比较\n');
fprintf('数值稳定性和计算效率的比较\n\n');

% 使用病态矩阵测试
epsilon = 1e-12;
A_ill = [1, 1, 1;
         epsilon, 0, 0;
         0, epsilon, 0;
         0, 0, epsilon];

fprintf('病态矩阵条件数: %.2e\n', cond(A_ill));

% Householder QR
[Q_house, R_house] = householder_qr_decomp(A_ill);
orthogonality_house = norm(Q_house' * Q_house - eye(size(A_ill, 2)), 'fro');

% 改进Gram-Schmidt (从前面的函数)
[Q_mgs, R_mgs] = modified_gram_schmidt(A_ill);
orthogonality_mgs = norm(Q_mgs' * Q_mgs - eye(size(A_ill, 2)), 'fro');

fprintf('正交性误差比较:\n');
fprintf('Householder: %e\n', orthogonality_house);
fprintf('改进GS: %e\n', orthogonality_mgs);
fprintf('Householder的优势: %.1f倍\n', orthogonality_mgs / orthogonality_house);

%% 6. 性能测试
fprintf('\n6. 性能测试\n');
fprintf('不同规模矩阵的计算时间比较\n\n');

sizes = [100, 200, 400, 800];
fprintf('矩阵规模    Householder    改进GS    MATLAB QR\n');
fprintf('--------    -----------    -------   ---------\n');

for i = 1:length(sizes)
    m = sizes(i);
    n = min(m, 100);
    A_test = randn(m, n);
    
    % Householder QR
    tic;
    [Q_h_test, R_h_test] = householder_qr_decomp(A_test);
    time_house = toc;
    
    % 改进Gram-Schmidt
    tic;
    [Q_mgs_test, R_mgs_test] = modified_gram_schmidt(A_test);
    time_mgs = toc;
    
    % MATLAB QR
    tic;
    [Q_matlab_test, R_matlab_test] = qr(A_test, 0);
    time_matlab = toc;
    
    fprintf('%6d      %9.4f s    %5.3f s   %7.4f s\n', ...
            m, time_house, time_mgs, time_matlab);
end

%% 7. 应用：最小二乘求解
fprintf('\n7. 应用：最小二乘求解\n');
fprintf('使用Householder QR分解求解最小二乘问题\n\n');

% 超定系统
m = 8; n = 3;
A_ls = [1, 1, 1;
        1, 2, 4;
        1, 3, 9;
        1, 4, 16;
        1, 5, 25;
        1, 6, 36;
        1, 7, 49;
        1, 8, 64];
b_ls = [2.1; 5.9; 13.1; 24.9; 41.2; 61.1; 84.8; 113.2];

fprintf('最小二乘问题: 多项式拟合\n');
fprintf('数据点数: %d, 多项式次数: %d\n', m, n-1);

% Householder QR求解
[Q_ls, R_ls] = householder_qr_decomp(A_ls);
x_house = R_ls \ (Q_ls' * b_ls);

fprintf('Householder QR解: [%.6f; %.6f; %.6f]\n', x_house);

% 计算残差
residual_house = A_ls * x_house - b_ls;
fprintf('残差范数: %.6f\n', norm(residual_house));

% 与MATLAB比较
x_matlab = A_ls \ b_ls;
fprintf('MATLAB解: [%.6f; %.6f; %.6f]\n', x_matlab);
fprintf('解的差异: %e\n', norm(x_house - x_matlab));

%% 8. 经济型QR分解
fprintf('\n8. 经济型QR分解\n');
fprintf('当m > n时，只计算前n列的Q矩阵\n\n');

A_econ = randn(10, 4);
fprintf('矩阵规模: %d×%d\n', size(A_econ, 1), size(A_econ, 2));

% 完整QR分解
[Q_full, R_full] = householder_qr_decomp(A_econ);
fprintf('完整QR: Q为%d×%d, R为%d×%d\n', size(Q_full), size(R_full));

% 经济型QR分解
[Q_econ, R_econ] = householder_qr_economic(A_econ);
fprintf('经济QR: Q为%d×%d, R为%d×%d\n', size(Q_econ), size(R_econ));

% 验证
fprintf('经济型分解误差: %e\n', norm(Q_econ * R_econ - A_econ, 'fro'));

%% 9. 更新和降级
fprintf('\n9. QR分解的更新\n');
fprintf('当矩阵增加一列时，如何高效更新QR分解\n\n');

A_base = randn(6, 3);
new_col = randn(6, 1);

fprintf('基础矩阵: %d×%d\n', size(A_base));
fprintf('新增列向量长度: %d\n', length(new_col));

% 原始QR分解
[Q_base, R_base] = householder_qr_decomp(A_base);

% 更新QR分解
[Q_updated, R_updated] = qr_update_column(Q_base, R_base, new_col);

% 验证
A_extended = [A_base, new_col];
[Q_direct, R_direct] = householder_qr_decomp(A_extended);

fprintf('更新方法与直接计算的差异:\n');
fprintf('Q矩阵差异: %e\n', norm(Q_updated - Q_direct, 'fro'));
fprintf('R矩阵差异: %e\n', norm(R_updated - R_direct, 'fro'));

end

function [v, beta] = householder_vector(x)
% HOUSEHOLDER_VECTOR 计算Householder向量
% 
% 输入: x - 列向量
% 输出: v - Householder向量, beta - 参数

n = length(x);
sigma = x(2:n)' * x(2:n);
v = [1; x(2:n)];

if sigma == 0
    beta = 0;
else
    mu = sqrt(x(1)^2 + sigma);
    if x(1) <= 0
        v(1) = x(1) - mu;
    else
        v(1) = -sigma / (x(1) + mu);
    end
    beta = 2 * v(1)^2 / (sigma + v(1)^2);
    v = v / v(1);
end

end

function [Q, R] = householder_qr_decomp(A)
% HOUSEHOLDER_QR_DECOMP Householder QR分解
% 
% 输入: A - m×n 矩阵
% 输出: Q - m×m 正交矩阵, R - m×n 上三角矩阵

[m, n] = size(A);
Q = eye(m);
R = A;

for k = 1:min(m-1, n)
    % 提取第k列的下半部分
    x = R(k:m, k);
    
    % 计算Householder向量
    [v, beta] = householder_vector(x);
    
    % 构造Householder矩阵并应用到R
    H_k = eye(m-k+1) - beta * (v * v');
    R(k:m, k:n) = H_k * R(k:m, k:n);
    
    % 更新Q矩阵
    Q(k:m, :) = H_k * Q(k:m, :);
end

Q = Q';  % 转置得到最终的Q

end

function [R, data] = householder_qr_compact(A)
% HOUSEHOLDER_QR_COMPACT 紧凑存储的Householder QR分解
% 
% 将Householder向量存储在R的下三角部分

[m, n] = size(A);
R = A;
beta_values = zeros(min(m-1, n), 1);
v_vectors = cell(min(m-1, n), 1);

for k = 1:min(m-1, n)
    x = R(k:m, k);
    [v, beta] = householder_vector(x);
    
    % 应用Householder变换
    H_k = eye(m-k+1) - beta * (v * v');
    R(k:m, k:n) = H_k * R(k:m, k:n);
    
    % 存储Householder数据
    beta_values(k) = beta;
    v_vectors{k} = v;
    
    % 在R的下三角部分存储v的后续元素
    if k < m
        R(k+1:m, k) = v(2:end);
    end
end

data.beta = beta_values;
data.v_vectors = v_vectors;

end

function Q = reconstruct_Q_from_compact(R, data)
% RECONSTRUCT_Q_FROM_COMPACT 从紧凑存储重构Q矩阵

[m, n] = size(R);
Q = eye(m);

for k = min(m-1, n):-1:1
    % 重构Householder向量
    v = [1; R(k+1:m, k)];
    beta = data.beta(k);
    
    % 应用Householder变换
    H_k = eye(m-k+1) - beta * (v * v');
    Q(k:m, :) = H_k * Q(k:m, :);
end

Q = Q';

end

function [Q, R] = householder_qr_economic(A)
% HOUSEHOLDER_QR_ECONOMIC 经济型Householder QR分解
% 
% 只计算前n列的Q矩阵

[m, n] = size(A);
Q = eye(m, n);
R = A;

for k = 1:min(m-1, n)
    x = R(k:m, k);
    [v, beta] = householder_vector(x);
    
    H_k = eye(m-k+1) - beta * (v * v');
    R(k:m, k:n) = H_k * R(k:m, k:n);
    Q(k:m, :) = H_k * Q(k:m, :);
end

Q = Q';
R = R(1:n, :);  % 只保留上三角部分

end

function [Q_new, R_new] = qr_update_column(Q, R, u)
% QR_UPDATE_COLUMN 向QR分解中添加一列
% 
% 输入: Q, R - 原始QR分解, u - 新列
% 输出: Q_new, R_new - 更新后的QR分解

[m, n] = size(R);

% 计算 Q^T * u
z = Q' * u;

% 使用Givens旋转消除z的下三角部分
R_new = [R, z];
Q_new = Q;

for k = n+1:-1:2
    if abs(R_new(k, n+1)) > eps
        % 构造Givens旋转
        [c, s] = givens_rotation(R_new(k-1, n+1), R_new(k, n+1));
        
        % 应用到R
        G = [c, s; -s, c];
        R_new([k-1, k], :) = G * R_new([k-1, k], :);
        
        % 应用到Q
        Q_new(:, [k-1, k]) = Q_new(:, [k-1, k]) * G';
    end
end

end

function [c, s] = givens_rotation(a, b)
% GIVENS_ROTATION 计算Givens旋转参数

if b == 0
    c = 1;
    s = 0;
else
    if abs(b) > abs(a)
        tau = -a / b;
        s = 1 / sqrt(1 + tau^2);
        c = s * tau;
    else
        tau = -b / a;
        c = 1 / sqrt(1 + tau^2);
        s = c * tau;
    end
end

end

function [Q, R] = modified_gram_schmidt(A)
% 改进的Gram-Schmidt (从前面复制)
[m, n] = size(A);
Q = A;
R = zeros(n, n);

for j = 1:n
    R(j, j) = norm(Q(:, j));
    if R(j, j) < eps
        error('矩阵A线性相关');
    end
    Q(:, j) = Q(:, j) / R(j, j);
    
    for k = j+1:n
        R(j, k) = Q(:, j)' * Q(:, k);
        Q(:, k) = Q(:, k) - R(j, k) * Q(:, j);
    end
end
end